{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Gamma regression for blood clotting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a preliminary example of this package's functionality, we provide an example of performing a Gamma regression, which is used when the response variable is continuous and positive. We have adapted the following canonical example of a Gamma regression from McCullagh & Nelder (1989). \n", "\n", "Nine different percentage concentrations with prothrombin-free plasma ($u$) and clotting was induced via two lots of thromboplastin. Previous researchers had fitted a hyperbolic model, using an inverse transformation of the data for both lots $1$ and $2$, but we will analyze both lots using the inverse link and Gamma family. \n", "\n", "The following initial plots hint at using a log scale for $u$ to achieve inverse linearity, as well as the fact that the two lots have different regression and intercept coefficients. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scikit_stan import GLM\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbsphinx": "hidden", "tags": [] }, "outputs": [], "source": [ "mpl.rc('axes.spines', top=True, bottom=True, left=True, right=True)\n", "#mpl.rc('axes', facecolor='white')\n", "mpl.rc(\"xtick\", bottom=True, labelbottom=True)\n", "mpl.rc(\"ytick\", left=True, labelleft=True)\n", "mpl.style.use('ggplot')\n", "\n", "\n", "# center images\n", "from IPython.core.display import HTML\n", "HTML(\"\"\"\n", "\n", "\"\"\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ATTRIBUTION: McCullagh & Nelder (1989), chapter 8.4.2 p 301-302\n", "bcdata_dict = {\n", " \"u\": np.array([5, 10, 15, 20, 30, 40, 60, 80, 100]),\n", " \"lot1\": np.array([118, 58, 42, 35, 27, 25, 21, 19, 18]),\n", " \"lot2\": np.array([69, 35, 26, 21, 18, 16, 13, 12, 12]),\n", "}\n", "bc_data_X = np.log(bcdata_dict[\"u\"])\n", "bc_data_lot1 = bcdata_dict[\"lot1\"]\n", "bc_data_lot2 = bcdata_dict[\"lot2\"]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l1, = plt.plot(bcdata_dict[\"u\"], bcdata_dict[\"lot1\"], \"o\", label=\"lot 1\")\n", "l2, = plt.plot(bcdata_dict[\"u\"], bcdata_dict[\"lot2\"], \"o\", label=\"lot 2\")\n", "\n", "plt.suptitle(\"Mean Clotting Times vs Plasma Concentration\")\n", "plt.xlabel('Normal Plasma Concentration')\n", "plt.ylabel('Blood Clotting Time')\n", "\n", "plt.legend(handles=[l1, l2])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l1, = plt.plot(bc_data_X, bc_data_lot1, \"o\", label=\"lot 1\")\n", "l2, = plt.plot(bc_data_X, bc_data_lot2, \"o\", label=\"lot 2\")\n", "\n", "plt.suptitle(\"Mean Clotting Times vs Plasma Concentration\")\n", "plt.xlabel('Normal Plasma Concentration')\n", "plt.ylabel('Blood Clotting Time')\n", "\n", "plt.legend(handles=[l1, l2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After this preliminary data analysis, we fit two lines to the two lots of data. Using $x = \\log u$, we fit a GLM to the data.\n", "\n", "The original results were as follows, and we recreate regression coefficients within a standard deviation of these values: \n", "\n", "$$\\text{lot 1:} \\quad \\hat{\\mu} ^{-1} = - 0.01655(\\pm 0.00086) + 0.01534(\\pm 0.00143)x $$\n", "$$\\text{lot 2:} \\quad \\hat{\\mu} ^{-1} = - 0.02391(\\pm 0.00038) + 0.02360(\\pm 0.00062)x $$\n", "\n", "As in previous work, we will fit two different linear models for each lot in the dataset. As usual, the $\\alpha$ parameter is the regression intercept and $\\mathbf{\\beta}$ is vector of regression coefficients and the parameter $\\sigma$ represents an auxiliary variable for the model. In this case, $\\sigma$ is the shape parameter for the Gamma distribution. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Initialize two different GLM objects, one for each lot. \n", "glm_gamma1 = GLM(family=\"gamma\", link=\"inverse\", seed=1234)\n", "glm_gamma2 = GLM(family=\"gamma\", link=\"inverse\", seed=1234)\n", "\n", "# Fit the model. Note that default priors are used without autoscaling, see the \n", "# API to see how to change these.\n", "glm_gamma1.fit(bc_data_X, bc_data_lot1, show_console=False)\n", "glm_gamma2.fit(bc_data_X, bc_data_lot2, show_console=False)\n", "\n", "print(glm_gamma1.alpha_, glm_gamma1.beta_)\n", "print(glm_gamma2.alpha_, glm_gamma2.beta_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As can be seen above, the fitted model has the following parameters, which are within one standard deviation of the results from past studies.\n", "\n", "$$\\text{lot 1:} \\quad \\hat{\\mu} ^{-1} = - 0.01437 + 0.01511 \\cdot x$$\n", "$$\\text{lot 2:} \\quad \\hat{\\mu} ^{-1} = - 0.02016 + 0.02301 \\cdot x$$\n", "\n", "As a verification of the accuracy of the fitted model, we can plot the fitted lines and the data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mu_inv1 = 1 /( glm_gamma1.alpha_ + glm_gamma1.beta_ * bc_data_X)\n", "mu_inv2 = 1 /( glm_gamma2.alpha_ + glm_gamma2.beta_ * bc_data_X)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mlot1, = plt.plot(bc_data_X, mu_inv1, \"r\", label=\"mu_inv lot 1\")\n", "mlot2, = plt.plot(bc_data_X, mu_inv2, \"b\", label=\"mu_inv lot 2\")\n", "l1, = plt.plot(bc_data_X, bc_data_lot1, \"o\", label=\"lot1\")\n", "l2, = plt.plot(bc_data_X, bc_data_lot2, \"o\", label=\"lot2\")\n", "\n", "plt.suptitle(\"Mean Clotting Times vs Plasma Concentration\")\n", "plt.xlabel('Normal Plasma Concentration')\n", "plt.ylabel('Blood Clotting Time')\n", "\n", "plt.legend(handles=[mlot1, mlot2, l1, l2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As this package is a wrapper around CmdStanPy, we can gather additional statistics about the fitted model with methods from that package. In particular, we can consider further statistics about the model by using CmdStanPy's summary method on the results of the fit. \n", "\n", "Notice that $\\mu$ (\"mu\") and the link-inverted $\\mu$ (\"mu unlinked\") are included as part of the model summary. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "glm_gamma1.fitted_samples_.summary()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "glm_gamma2.fitted_samples_.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additional information about the model and various visualizations can be revealed by Arviz, which seamlessly integrates with CmdStanPy components. Consider the following. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import arviz as az\n", "az.style.use(\"arviz-darkgrid\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "infdata = az.from_cmdstanpy(glm_gamma1.fitted_samples_)\n", "infdata" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "az.plot_posterior(infdata, var_names=['alpha', 'beta', 'sigma']);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "az.plot_trace(infdata, var_names=['alpha', 'beta', 'sigma'], compact=True);" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.8.5 64-bit", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" }, "vscode": { "interpreter": { "hash": "9650cb4e16cdd4a8e8e2d128bf38d875813998db22a3c986335f89e0cb4d7bb2" } } }, "nbformat": 4, "nbformat_minor": 4 }